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The melting transition of a system of hard cubes is studied numerically both in the case of freely 
rotating cubes and when there is a fixed orientation of the particles —parallel cubes. It is shown that 
freelly rotating cubes melt through a first-order transition, whereas parallel cubes have a continuous 
transition in which positional order is lost but bond-orientational order remains finite. This is 
interpreted in terms of a defect-mediated theory of melting 



I. INTRODUCTION 

Freezing is probably the most unavoidable phase tran- 
sition of a classical system of identical particles, occurring 
when temperature is reduced sufficiently. The most im- 
portant difference between the low temperature (solid) 
phase and the high temperature (fluid) phase is that in 
the solid phase long range correlations between the coor- 
dinates of particles exist, whereas in the fluid phase this 
correlation is short range, decaying exponentially with 
distance. Freezing (or melting) is —for three dimensional 
systems— a first-order transition, with discontinuities in 
the first derivatives of thermodynamic potentials at the 
transition point. Compared to the liquid-gas transition 
or magnetic transitions in model systems, the advance in 
the theoretical understanding of the melting transition 
has been considerably slow. 

Historically, there have been two main theoretical 
approaches, to the problem of melting. Pcrturbative 
arguments^ starting from the solid phase are able to pre- 
dict that the solid structure will be unstable when tem- 
perature is increased sufficiently, but cannot automati- 
cally predict the characteristics of the fluid phase. De- 
scriptions starting from the fluid phase (virial expansions, 
Ornstcin-Zernicke equations, etc.), although give very ac- 
curate descriptions of the fluid phase in some model sys- 
tems, usually do not predict at all a transition to a solid 
phase. The most successful account of the melting transi- 
tion arising from this-line of thinking is an order parame- 
ter theory of meltingcl in which the most stable structure 
at a given temperature is obtained by minimizing a free 
energy functional which depends on the intensity of the 
Bragg peaks of the structure. These intensities are zero 
in the fluid phase, and finite in the solid phase. 

A more recent approach to the problem is the theory 
of defect-mediated melting transition originally proposed 
for two-dimensional (2D) systems. The original idea is 
the following.uu The melting from the solid to the fluid 
phase is a consequence of the proliferation of defects in 
the perfect crystalline structure that exists at zero tem- 
perature. These defects, called dislocations, appear in 
pairs defect-antidefect, and have a finite binding energy. 



As long as the size of the pair is finite, the positional 
order of the system is only perturbed at small length 
scales. When temperature is increased beyond a critical 
value T m , the size of the pair becomes infinite, or in other 
words, free dislocations can exist in thermal equilibrium 
in the system. The existence of free dislocations disrupts 
the long range positional order of the system, and drives 
the melting transition of the crystal. The transition is. 
predicted to be continuous, of Kosterlitz-Thouless type.cl 
It was further realized that a crystal with a finite con- 
centration of free dislocations does not behave exactly as 
a usual fluid. In fact, dislocations destroy the positional 
order of crystals, but not the orientational order that can 
still be present. In this scenario, the state of the system 
for T > T m corresponds to a fluid .with orientational or- 
der, known as the hexatic phase.B Within the hexatic 
phase a new kind of defects appear that destroy the ori- 
entational order at high enough temperatures. These de- 
fects were called disclinations by Halperin and Nelson. 
At a certain temperature T Q disclination pairs unbind, 
and the system transforms in a usual fluid through a new 
Kosterlitz-Thouless transition. The critical values T a and 
T m depend on the self energies and interaction energies of 
dislocations and disclinations, which in turn are depen- 
dent on the parameters of the model. It should be noted 
that a necessary condition for this two-step melting pro- 
cess is that T Q > T m , since orientational order can exist in 
the absence of positional order, but positional order can- 
not exist in the absence of orientational order. If the cou- 
pling between dislocations and disclinations is strong, the 
two continuous transitions can merge into a single first- 
order transition. trQ Many theoretical and experimental 
work was since then devoted to the confirmation of this 
theory of two-dimensional melting (that is known as the 
Kosterlitz-Thouless-Halperin-Nelson- Young — KTHNY— 
theory) and it was found in fact that the melting transi- 
tion in two dimensions may be. a single first order ,B or a 
two-step continuous transitional depending on the system 
studied. 

Due to its usefulness in 2D systems, it is tempting 
to apply the defect-mediated theory of melting to three- 
dimensional (3D) systems!] In three dimensions, disloca- 
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tions are one dimensional defects that form closed loops, 
or open lines that begin and end at the surfaces of the 
sample. In a 3D solid at low temperatures only small 
dislocation loops exist, but beyond a critical temperature 
infinitely large loops are present at equilibrium and the 
solid looses its positional long-range order. This picture 
of the transition is closely related to the superconducting- 
normal transition driven Jay the proliferation of vortex 
loops in superconductors ,L2I or the superfluid- normal liq- 
uid transition in He 4 Jill Again, the transition driven by 
dislocations alone is continuous and the fluid above melt- 
ing still has orientational long range order that generates 
a residual resistance to torsion, not present in a normal 
fluid.lla In analogy to the situation for 2D systems, a new 
class of defects responsible for the disappearance of ori- 
entational order must be introduced, and a possible new 
transition at higher temperatures should be expected. 

From this point of view, the striking feature of the 
melting of identical particles in three dimensions is 
that it is always a first-order transition. This may in- 
dicate that the disclination-unbinding and dislocation- 
unbinding transitions in three dimensions are strongly 
coupled, in such a way that they promote each other and 
make the transition first-order, but why this is always so 
is not known. The existence of a model system that melts 
through a continuous transition into a fluid with orien- 
tational order would be of importance in giving insight 
into the KTHNY theory in three dimensions. 

The aim of this work is to analyze numerically a simple 
model which displays a continuous melting transition in 
three dimensions. The model is a system of impenetra- 
ble cubes, which have fixed orientation in-space, the same 
for all cubes (parallel hard cubes, PHC).E3 KirkpatrickEil 
showed that this model has a continuous transition to a 
simple hypercubic solid structure in infinite dimensions, 
and suggested that this would also be so in three dimen- 
sions. The two main reasons to expect a continuous melt- 
ing for PHC are the following. First, a cubic crystal lacks 
—in a Landau description of its melting— a third order 
term in the free energy functional that would favor the 
transition to be first-order .113 This kind of terms appear 
for crystalline structures that possess three Bragg vectors 
Gi, G2, G3 lying on the first maximum of the diffraction 
pattern, and satisfying the relation Gi + G2 + G3 = 0. 
These vectors do not exist for a simple cubic structure. In 
addition, bond-orientational ordered will be strongly en- 
hanced in PHC compared for instance to spheres because 
a fixed orientation of each cube favors a neighborhood 
in which cubes arrange with the same orientation. This 
raises the possibility for the orientational order to persist 
up to higher temperatures than the translational order. 
For comparison, the case of freely rotating hard cubes 
(FRHC) will also be studied, and it will be shown that 
in this case the melting is a usual first-order transition 
into an isotropic fluid. 



II. NUMERICAL TECHNIQUE AND RESULTS 

The numerical method used to simulate the system is a 
standard Monte Carlo-Metropolis algorithm in the NPT- 
ensemble. The positions of the cubes are characterized 
by the coordinates of their centers. A trial movement 
of a particle consists of a displacement to a new posi- 
tion chosen randomly inside a cube of linear size 0.01 I 
centered at the old position (/ is the linear size of the 
particles). The new position of the particle is accepted 
as long as there is no overlap with any other particle. 
After all particle coordinates are updated a trial global 
rescaling of all particle coordinates and system size by 
a factor within the range 1 ± .01 is proposed. If this 
change does not produce particle overlapping, then it is 
accepted according to the Metropolis algorithm with an 
energy change dE given by dE = PAV - Nk B TAV/V 
(N is the total number of cubes, P is pressure and V 
is the volume of the system). In the case of FRHC, in 
addition to the center-of-particle coordinates, the three 
Euler angles are necessary to characterize the position 
of each cube. These angles are updated at the same 
time as coordinates, the elemental change in each step 
is chosen to be ~ 0.1. Since there is no configurational 
contribution to the energy of the system, the equation 
of state depends only on the relation T/P. All results 
are presented as function of the adimensional tempera- 
ture T* = kBVQ 1 (T/P) (that will be refered to simply as 
"the temperature"), where vq = I is the volume of each 
cube.li3 

The zero temperature state of the system of cubes 
(both parallel or freely rotating) is highly degenerate, be- 
cause along any of the main crystalline directions, rows 
of cubes can be displaced an arbitrary amount without 
changing the volume of the system. However, at finite 
temperatures the cubic configuration with long range po- 
sitional order has larger entropy than any row-displaced 
configuration, and thalhermodynamically stable state is 
a simple cubic lattice .113 Even for the small systems that 
we are going to simulate, this entropy is greater than 
the one that can be gained by displacing rows of cubes 
(which is of the order of ln(iV) /N), and configurations 
with displaced rows never show up in the simulations in 
the temperature ranges of interest. 

When temperature is increased sufficiently the crystal 
melts. This melting is qualitatively different for PHC and 
for FRHC. In the case of FRHC the melting occurs via 
a standard first order transition. Results of simulations 
are presented for a system of 125 particles. The system 
was initialized in a perfect cubic structure at low temper- 
ature, and a simulation was performed by increasing and 
then decreasing temperature. At each temperature 5000 
Monte Carlo steps were used for thcrmalization and then 
20000 steps were used to compute quantities of interest. 
In Fig. 0(a) we see the evolution of the inverse packing 
fraction v = V/ (Nvq) of the system. It shows a clear hys- 
teretic behavior indicating a first-order phase transition 
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for T* ~ 0.15, where density changes between ~ 0.45 
and ~ 0.52. Also shown in this figure as a dotted lins 
are the values predicted from a cell-theory for the solid,EJ 
which gives a reasonable approximation to the real equa- 
tion of state up to the melting temperature. It would be 
nice to have the expressions for the virial coefficients of 
FRHC to fit the fluid part of the curve, but these are not 
available to the required order to get a good fitting. For 
comparison, in Fig.fl .the Carnahan-Starling equation of 
state of hard spheres!] is shown. The only free parameter 
is the sphere volume that was chosen to be 1.2w - 




FIG. 1. Inverse packing fraction v (a), Bragg intensities B c 
(b, full symbols) and nearest neighbors orientational order B2 
(b, open symbols) as a function of the adimensional tempera- 
ture T* for a system of 5 x 5 x 5 FRHC, upon cooling (starting 
from a disordered configuration) and heating (see text for def- 
inition, £>2 and B c are given in arbitrary units). In (a) the 
dotted line is the prediction from a cell theory of the solid, 
and the dashed line is the behavior of hard spheres with an 
effective volume of 1.2i>o in the fluid phase. 

In Figjl](b) two different indicators of the order in the 
system support the conclusion that the melting transition 
of FRHC is first order. The parameter B c is extracted 
from the diffraction pattern of the structure, and it is 
defined as 

4 f 

B c = I / D (M, <p) Y*,m (0, <p) x 

x 5 (k -h)k 2 dk sin (0) d0dip\ 2 , (1) 

where D (fc, 6, (p) is the intensity of the diffraction pat- 
tern in polar coordinates, the delta factor picks up the 
values at the first maximum of the diffraction pattern 
(k% = 27TV 1 / 3 ), and the spherical harmonics Y^ m collect 



the part with cubic symmetry of the diffraction pattern. 
The value of B c is different from zero if the system pos- 
sesses long range positional order S3 

The relative ordering of neighbor particles B2 is defined 

as 



B 2 = £ 



D 2 (r, 9, <p) K (r) Y 4>m (fi, <p) x 

xr 2 dr sin (9) d9dtp \ 2 (2) 



with D 2 (r, 9, if) being the pair distribution function of 
particles at distance r, along the spatial direction (9, ip) . 
The kernel K (r) cuts off the integral beyond some dis- 
tance. The results are qualitatively insensitive to the 
exact form of K (r) , in the results presented below K (r) 
was taken to be 1 for r < 1.5U 1 / 3 , and for r > 1.5U 1 / 3 . 
The value of B 2 is different from zero if the system pos- 
sesses long range orientational order. 

All these indicators of ordering vanish at the melting 
transition, with the same hysteretic behavior as that of 
the volume. The unambiguous determination of a first- 
order phase transition would require the study of the 
volume histogram at the transition temperature, which 
should have a double peak structure associated to the co- 
existence of a solid and a fluid phase. Unfortunately, the 
simulation of FRHC is very time consuming so as to carry 
out this program. Partial checks were performed, how- 
ever. In a simulation around the transition temperature 
(T* = .15) the volume of the system stabilized around 
different values, depending if the initial configuration of 
the system was chosen random or ordered. These values 
were the ones expected from Fig. fl](a). Partial simu- 
lations in systems up to 512 particles were performed, 
and the results are consistent with a first-order melting 
transition for FRHC. 

If the cubes are restricted to be parallel to each other, 
the nature of the melting transition changes qualitatively. 
Results of simulations for this case are shown in Fig. |^ 
for a system of 216 particles. The volume of the sys- 
tem does not show any abrupt change, but a continuous 
and reversible (on heating and cooling) behavior. The 
parameter characterizing the crystalline order B c dimin- 
ishes strongly around T* — 0.4, where the system has a 
density ~ 0.5, suggesting a continuous melting. The lo- 
cal orientational order (characterized by B 2 ), in spite of 
decreasing near the transition remains finite at high tem- 
peratures. This characteristic is not surprising since the 
orientational order is favored by the equal orientation of 
all cubes. In Fig. |we can see also the predictions for the 
volume from the lowest order cell model of the solid and 
the seventh order virial expansion for the fluid.Ej These 
expressions give a good approximation to the simulated 
values for all temperatures. 
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FIG. 2. Same as Fig.[j] for a system of 6 x 6 x 6 PHC. In 
(a) the dotted line is the prediction from a cell theory of the 
solid, and the dashed line is the equation of state for PHC to 
seventh order virial expansion. 

If the melting of PHC is really a second order phase 
transition, the behavior of the order parameter of the 
transition (that can be taken to be the crystalline order 
parameter B c ) must obey scaling laws as a function of 
the system size. In particular, different simulations of B c 
in systems of different sizes-,!/ (= V 1 ^ 3 /!) must obey a 
scaling relation of the formEil 



7 (V 



T*)L l ' v 



(3) 



where / is a universal function, v and \x are two critical 
exponents and is the thermodynamical melting tem- 
perature. The exponent v characterizes the divergence at 
the thermodynamic melting temperature of the cor- 
relation length. The result of simulations for systems of 
216, 512, and 1000 particles are shown in Fig. || The vol- 
ume and the orientational order Bi show no detectable 
dependence on size, whereas the crystalline order B c has 
a clear size dependence. Results for B c for different sys- 
tem sizes can be collapsed reasonably well onto a single 
curve when plotted as B C L^ vs (T* - T£) L 1 /", with pa- 
rameters T* t = 0.40 ±0.02, /i = 4.0 ±0.5, v = 0.50 ±0.05. 
This value of v is lower than the one corresponding to a 
three dimensional XY model, or the loop model for the 
normal-to-superconducting transition {y — .666 ± .003) 
that is supposed to be in the same universality class of 
our model if the melting can be described by the KTHNY 
theory. However, to be able to unambiguously decide 
this point, more simulations in larger systems are needed. 
The density of the system at melting is 0.48 ± 0.02. 
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FIG. 3. Inverse packing fraction v (a), orientational order 
B2 (b), and crystalline order B c (c) as a function of tempera- 
ture for systems of PHC of different sizes. The values shown 
correspond to an average upon cooling and heating. In (d) 
the curves of B c are scaled according to a second order phase 
transition using the adimensional linear size of the system 
L = V 1/3 /l. 

The fluid formed by the parallel cubes above melting 
is not a usual isotropic fluid. This is obvious since some 
spatial orientations are singled out by the particular form 
of the particles. The fluid phase of PHC is the analogous 
of the hexatic phase of the KTHNY theory. Within this 
framework, the difference between FRHC and PHC is 
clear: parallel cubes keep the long range orientational 
order even when positional order has been lost, and the 
KTHNY theory predicts a continuous melting if only po- 
sitional order is lost. FRHC have the possibility of loos- 
ing both positional and orientational order, and this is in 
fact what happens at a unique temperature in a discon- 
tinuous form. 

It may be of interest to compare the_fluid of PHC 
with the nematic phase of liquid crystals Ej In that case, 
molecules orient along a preferred spatial direction (i.e., 
they possess molecular-orientational order). Upon cool- 
ing, this structure transforms usually into a smectic-A 
phase in which a long range positional order is established 
along the direction characterizing the nematic phase. 
This transition may be first or second order depending 
on the material. At a lower temperature the smectic- 
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A phase can undergo a transition to a crystalline phase. 
In our case, parallel cubes single out three orthogonal 
and equivalent directions in space, and upon cooling the 
system freezes into a solid phase, with crystalline order 
in all directions. The melting of PHC has no analogous 
in the transitions that occur in liquid crystal systems. 
Note that for the case of cubes the oriented phase has to 
be stabilized from outside, whereas the nematic phase in 
liquid crystals may be generated by molecular hard core 
interactions only. 



III. CONCLUSION 

In summary, I have shown numerical results on a sim- 
ple model that displays a continuous melting transition 
in three dimension —namely a system of parallel hard 
cubes. The melting of this system can be qualitatively 
interpreted in terms of the KTHNY theory of defect- 
mediated melting. The melting temperature was esti- 
mated to be = 0.40 ± 0.02, and the critical density 
is 0.48 ± 0.02. The critical exponent of the correlation 
length is v = 0.50 ± 0.05. At the melting transition only 
positional order is lost, orientational order remains finite 
because it is favored by the geometric form of the par- 
ticles. If the cubes are allowed to rotate, the melting is 
a usual first order transition where both positional and 
orientational order are lost. 
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